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Abstract 


■t- 


In  this  project,  the  main  objective  was  to  develop  multiresolution  wavelet  algo¬ 
rithm  to  study  the  flame  acceleration  and  to  understand  the  mechanism  of  the  tran¬ 
sition  of  deflagration  to  detonation.  In  the  first  half  year  of  the  project,  The  author 
completed  the  study  of  two  dimensionail  detonation  waves  based  on  hybrid:  high  order 
methods  (a  combination  of  high  order  Essentially  Non-oscillatory  (ENO)  methods  and 
spectral  methods  and  Shock  Tracking  methods).  The  result  on  the  detonation  waves 
was  published  in  the  AIAA  Journal.  In  the  remaining  time  of  the.two  iye^s,  first, 
the  author  completed  the  theoretical  and  algorithmic  studies  of  the  adaptive  wavelet 
method,  which  could  handle  nonperodic  boundary  conditions  and  nonlinear  time  de¬ 
pendent  PDE’s  (The  result  was  published  in  the  SIAM  Journal  of  Nuinerical  Analysis). 
Secondly,  the  author  implemented  the  adaptive  wavelet  methods  for  the  solution  of  one¬ 
dimensional  flame  propagation;  thirdly,  the  author  developed  a  Fortran  code  WL2D 
(more  than  13,000  lines)  for  the  two  dimensional  multi-scale  wavelet  algorithms  with 
an  efficient  data  structure  and  implemented  a  second  order  Implicit  fi^torized  scheme 
for  the  adaptive  wavelet  methods. 
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1  Hybrid  High  Order  IVlethods  for  iDetohation  Waves 

Partially  supported  by  this  grant,  we  completed  the  development  of  high 'order  hybrid  nu¬ 
merical  simulation  of  two  dimensional  detonation  waves.  The  major  finding pf|  this  work  was 
that  the  cellular  structure  of  detonation  waves  depended  very  sehsitiyely  ,dn  the  numerical 
dissipation  of  the  algorithms  representing  detonation  fronts.  Further  studies  of  the  work  in 
three  dimensional  cases  were  needed  to  understand  the  three  dimensional  effects  of  detona¬ 
tion  waves.  One  paper  summarizing  the  results  of  this  study  was  published  in  AIAA  Journal, 
Vol.  33,  Number  3,  pp  1248-1255. 

2  Parallel  Multi-scale  Wavelet  Algorithms 

In  an  attempt  to  design  multiscale  methods  for  the  study  of  deflagratibn'  tojdetonation  tran¬ 
sition  (DDT)  problem,  we  constructed  a  wavelet  collocation  multi-resoluiibn  algorithm  for 
the  initial  value  boundary  problem  of  nonlinear  PDF’s.  The  key  obmppiieht  in  this  colloca¬ 
tion  method  was  a  so-called  “Discrete  Wavelet  Transform”  (DWT)  whibh  niaped  a  solution 
between  the  physical  space  and  the  wavelet  coefficient  space.  The  DWT,  transformation  only 
took  O(NlogN)  operations  where  N  was  the  total  number  of  unknowns.  Thbrefore,  the  non¬ 
linear  term  in  the  PDE  could  be  easily  treated  in  the  physical  spape,  and  ihe.  derivatives  of 
those  nonlinear  terms  then  computed  in  the  wavelet  spane.  The  wavdet  collpcation  methods 
had  the  following  advantages:  (a)  the  capability  to  handle  arbitra.iy;jiiph.-periodic  boundary 
conditions;  (b)  the  capability  to  treat  general  nonlinearity  through,  the  cbllocations  of  the 
PDF’s;  (c)  flexibility  of  adaptive  meshing  in  regions  where  high  gradients  solution  occur 
as  in  shock  waves  and  turbulent  premixed  flames;  (c)  Capability  pf  paraliei;  and  multi-grid 
implementations  . 

The  following  paragraph  summarizes  the  key  technical  details  of  the- Algorithms.  One 
paper  containing  some  of  the  details  of  the  results  of  this  study  was  phbiished  in  the  Journal 
of  SIAM  Numerical  Analysis,  Volume  33,  Number  3,  pp.  937*970,  June:  19.96. 

2.1  Discrete  Wavelet  Transform  (DWT) 

The  Discrete  Wavelet  Transformation  (DWT)  maps  between  a  function  at  its  sample  points 
to  the  coefficients  of  its  wavelet  interpolation  expansion.  This  transfprrn  in  both  direction 
will  only  take  0{NlogN),  N  =  -f  L  -  1  operations  for  the  !ff^(/)7wiayelet  basis. 

Let  f{x)  in  Hq{I),  and  we  intend  to  construct  a  wavelet  interpolatipn 

Vj/ix)  e  Vo  ®  Wo®  for  J  >  b-  S.; 

First  introduce  two  index  spaces,  for  Wj,j  >  0  we  define 

/Cj  =  {-l,0,---,nj -2}  "iV:  (1) 

where  Uj  =  DirriWj  =  2^L  and  for  Vq  we  define  ..  -i.^.  : 

x:_:  =  {-i,o,-v,T~3}  ■  (2) 
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P2 


Second,  define  the  collocation  points  in  domain  I  =  [0,  L]  associated  with  wavelet  space 
>  0  by 


(j)  _  k  +  1.5 


= 


2^ 


(3) 


and  for  Vo  the  collocation  points  are 

^  +  2,  k  ^  K-i  (4) 

The  interpolation  operator  Ivof  in  Vq  interpolates  f(x)  on  collocation  points 
while  the  interpolation  operator  in  Wj  interpolates  a,ny  functiph^  j^ir)  in  on 

collocation  points  The  construction  of  the  interpolant  involves  the 

inversion  of  a  tridiagonal  matrix  with  size  dimWj  =  nj  ot  dimVh  =  X  —  ly  Tesp 

Now  let  us  assume  that  the  values  of  a  function  /(x)  €  ffoU)  ire?. given  on  all  the 
collocation  points  6  <  j  <  J,  then  the  wavelet  interpolation  Vjf{x)  G 

Vo  0  IVo  ©  Wi  “  •  •  ©  Wj  for  J  >  0  is  defined  as 


such  that 


where 


and  for  j  >0, 


‘Pj/iXt’)  =  f(Xk^),  for  k  e  <  j 

/-i(l)=Iv./(i)=  E  /-i,i«0.tW  €V4  ,  V? 
Jk€/C_i 


(5) 


(6) 


fj{x)  =  =  x:  €  WO:  . :  . 

Let  us  denote  f  =  the  values  of  f{x)  on  ail  interpolation  points  on 

all  levels,  =  {/(xjf^)}igCj-,  j  >  —1  and  f  =  the’Wa,velet  coefficients 

in  the  expansion  (5),  > -1. 

A  fast  discrete  wavelet  transform  (DWT)  was  proposed  for  an,  el&cient  transformation 
between  point  values  f  and  f  with  operation  counts  of  order  0(iVhifif)iin  both  directions  . 
The  efficiency  of  the  transformation  is  based  on  the  fact  that  fc  €  Kj  form 

a  hierarchical  nodal  basis  on  all  levels  of  collocation  points  {x*  —1  <  j,  <  J,  k  £  Kj. 

2.2  Wavelet  Collocation  Methods  for  PDB’s 

We  consider  a  collocation  method  based  on  the  DWT  transform  for  tinie  dependent  PDE’s. 
Let  u  =  u(x,  t)  be  the  solution  of  the  following  initial  boundary  value  problem 

Ut  +  fx{u)  =  Uxx  +  ff(u),  X  €  [0,  i],.t  >  0  .  ,.r,  ,w  '■  r 
u(0,t)  =5o(t)  ■  ■  ■• 

Bu{L,t)  =5i(t)  '  •  ■ 

Bu(i,0)  =/(x) 


(7) 
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where  B  is  the  boundary  condition  operator  which  could  be  either  Dirichlet  boundary  con¬ 
ditions  or  Neuman  type  or  Robin  type  boundary  conditions. 

The  numerical  solution  uj{x,i)  will  be  represented  by  a  uni (jue  decomposition  in  the 
cubic  spline  wavelet  decomposition  of  -  Vo  ®  ®  •  0  Wj^  J  >  Oj  namely 

uj(a:)  =  Ij,ju(i) -f  7/_i(x)-|-uo(x)  +  -f  uj(x)  (8) 

where  cubic  spline  consists  of  the  nonhomogenuity  of  u(i,  f)  on  both  boundaries, 

and  the  coefficients  axe  all  functions  of  t.  Using  the  DWT  transform,  we  can  also 

identify  the  numerical  solution  uj{x,t)  by  its  point  values  on  all  collocation  points,  we  put 
all  these  values  in  vector  u  =  u(t),  i.e. 

u  =  u(t)  =  -  •  • , 

where  € /Cj- for  j  >  — 1. 

To  solve  for  the  unknown  solution  vector  u(t),  we  collocate :the  PDE  (7)  on  all  collocation 
points,  then  we  have  the  following  semi-discretized  wavelet  collocation  niethod. 

Semi-Discretized  Wavelet  Collocation  Methods 

'  fx{uj)  = 

Buj{0,t)  =go{t) 

I  Buj{L,t)  -gi{t) 

[  uj{x]^\0)  =/(xSf^) 

where  fc  €  Xlj  —  1  <  j  <  J. 

Computation  of  /z(xt^)  =  fxiuj{x]l^) 

Step  1  Given  u  =  (u(“^\uW,  •  •  -  ^  compute  fc  6  /Cj,jf  > 

define 

f  =  (f(-l)^fW 

Step  2  Compute  the  wavelet  interpolation  expansion  using  DWT  transform  for  f  as  in  (5) 

Step  3  Differentiate  the  interpolation  expansion  and  evaluate  at  all  collocation  points  x]i\ 

which  is  taken  as  fzi(xk^))  =  fx{'Uj{xk^)). 

The  total  cost  of  computing  the  derivatives  will  be  (5J  -|-  12)A^  <  ;5Nl6gN. 

2.3  Adaptive  Meshing 

In  equations  (8),  uj{x)  is  expressed  using  the  full  set  of  collocation  points  {x^-'^}.  As  most 
of  the  wavelet  expansion  coefficients  uy,*  for  large  j  can  be  ignored  within  'a  given  tolerance 
£.  So  we  can  dynamically  adjust  the  number  and  locations  of  the  cpllociatibn  points  used  in 
the  wavelet  expansions,  reducing  significantly  the  cost  of  the  scheme  vsrhile  providing  enough 
resolution  in  the  regions  where  the  solution  varies  significantly. 
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(9) 


— 1  and 


Let  e  >  0  be  a  prescribed  tolerance  and  j  >  0,  £  =  £(€).  =  —  log  e/  log  a),  a  = 

13.98.  '  • 

Step  1.  First  locate  the  range  for  the  index  k,  ■ 

=  (10) 

such  that 

i  =  (11) 

Step  2.  Ignore  uj^k  in  (8)  for  fci  <  A:  <  /.-, ?  =  1,  •  -  • ,  m, k{  =  k-  4-  +  3, /<  =  ^  -  3, 

namely  we  redefine  Uj(i)  as 

uj(x):=  ^ 

igXTyXAC; 

where  ;C'  = 

Step  3.  The  new  collocation  points  and  unknowns  will  be 

€  x:_,  if  J  =  -i;k  €  /Cj\/c;  i£j  >’o. 

2.4  Data  Structure  for  Two  Dimensional  Adaptive  Wavelet  Algo¬ 
rithms 

We  have  developed  an  efhcient  data  structure  for  two  and  three  dimensional  wavelet  ap¬ 
proximations.  Being  an  adaptive  scheme  in  nature,  an  efficient  data  stihcture  for  handling 
numerical  solutions  in  a  wavelet  framework  is  the  first  step  in  m^ihg  a  theoretical  poten¬ 
tial  into  real  application.  We  have  used  speirse  matrix  data  structure  (cpihpressed  row  and 
compressed  column  vector  technique)  in  treating  the  data  structure  oii  e^  wavelet  space 
Wf  X  WJ  in  the  case  of  two  dimensional  approximations.  Such  approach  have  the  advantage 
of  only  storing  the  mesh  points  used  in  the  adaptive  wavelet  approximations  and  easy  access 
to  numerical  data  on  each  constant  x  and  y  lines  for  the  numericail  differentiations. 

Data  Structure  -  1-D  . 

One  Dimensional  Case: 

u  =  •  •  • ,  • '  1 

€  V'„©(ro0w,0.-iyj 

where  rij  -  number  of  mesh  points  on  level  Wj 

[1]  pointer{j)  ~  pointer  of  first  element  of 

[2]  npts(j)  -  number  of  elements  in  —  rij 

[3]  mdea:(l  :  rij)  -  collocation  point  location  indices  ^  J  ' 
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Data  Structure  -  2-D 


X  ,  ly?'  =  •••  W?  X  ly?  ••• 

^;x=5-l  jz  ^  jz  ^  jy 

Mesh  on  level  {jx,jy)  —  x 

•  Compressed  Row  Form  ; 
i*th  row: 

lenrow{i,jx,jy)  -  length  of  i-th  row 
ipr{i,jxjy)  -  pointer  of  1st  mesh  on  i-th  row 
icri(:,  -  column  indices  of  i-th  row 

•  Compressed  column  Form  : 
i-th  column; 

lencol{i,jx,jy)  •  length  of  i-th  row 
ipc{i,jx,jy)  -  pointer  of  1st  mesh  on  i-th  row 

irn{:,jx,jy)  -  column  indices  of  i-th  row  v  ^ 

2.5  Fast  Time  Integration  -  Factorized  ADI  Wavelet  Approach 

On  the  issue  of  time  integration,  we  have  successfully  iinplemented  a  second  order  implicit 
factorized  scheme  of  Beam  and  Warming  type  for  the  adaptive  wavelet  ttiethods.  Unlike  many 
other  adaptive  methods  (finite  element  and  finite  difference),  the  idaptivei-wavelet  methods 
actually  can  be  easily  implemented  using  an  ADI  approach..?  So,i  solution':  of  the  algebraic 
systems  from  the  implicit  discretization  of  'Z-dimensional  diffusioii  operators  is  replaced  by 
that  of  only  1-dimensional  operators,  thus  speeding  up  the  time  integfatmUs  tremendously. 
Example:  Beam  Warming  schemes  for  Euler  Equations. 

du  dF^  dG^_.  ■■  '  . 

dt^dx^dy~ 

(/  +  ^D,A’'){I  +  y  i),B")u"+'  , 

=  (/  - 

D®  and  the  wavelet  derivatives. 
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3  Computational  Results  of  Multi-scale  Wavelet  Algo¬ 
rithms  for  2-D  flames 

We  have  developed  a  two  dimensional  code  WL2D  based  on  the  adajitive:  wavelet  approx¬ 
imations.  The  following  eight  pages  summarize  the  simulation  results  of  two  dimensional 
cellular  and  planner  flame  propagations. 

[1]  Two  Dimensional  Perturbed  Cellular  Flame 


0t  =  A0  — K— -|-fl 
r 

AC  Cr 


Ci  =  — 

L  T 

where  A  is  the  Laplace  Operator  Initial  Condition: 


1 


C  =  (l-G)  +  0(^) 


-|-  perturbation  of  flame  fronts 
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[2]  Two  dimensional  Perturbed  planar  flames 


0,. 

Ct 


=  A0  +  n 


where  A  is  the  Laplace  Operator  Initial  Condition: 


0  =  exp  z  if  I  <  0 
=  1  if  z  >  0 
C  =  1-0 


■ 

'  '  ‘‘v' 
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Reactant  Mass  Fraction 
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